//================================================================================================================================//
#include "mpi.h"
#include "silo.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "ctype.h"

#define n_var_total 33
#define char_len_5 5
#define char_len_128 128
//================================================================================================================================//

int main(int argc, char **argv)
{
//--------------------------------------------------------------------------------------------------------------------------------//
// ro, gp, mp, te, beta, vx, vy, vz, vv, bx, by, bz, bv, jx, jy, jz, jv, fgx, fgy, fgz, fgv, fggx, fggy, fggz, fggv, 
// fmgx, fmgy, fmgz, fmgv, fmtx, fmty, fmtz, fmtv
//--------------------------------------------------------------------------------------------------------------------------------//
   FILE *datamap_unit;
   DBfile *dbfile_unit;
   DBoptlist *opt_list;
   DBmrgtree *mrgtree;

   char i_t_char[char_len_5], i_cpu_char[char_len_5], block_id_char[char_len_128];
   char directory[char_len_128], datavis_name[char_len_128], char_tmp[char_len_128];
   char **var_name, **coord_name, **vector_comp_name, ***scalar_name, ***vector_name;
   char **mesh_name, *mesh_name_block;
   char *level_name[1];

   int db_var_type;
   int i_var, j_var, n_var, len_str, i_cpu, i_level;
   int i_grid, i_grid_tmp, i, j, k, i_tmp, j_tmp, k_tmp;
   int i_e, j_e, k_e;
   int p1, p2, p3, p4, p5, p6, p7;
   int nx_no_bnd, nx_block, nx_total;
   int ny_no_bnd, ny_block, ny_total;
   int nz_no_bnd, nz_block, nz_total;
   int nxny, nxnynz, nxnyb, nxnynzb;
   int ini_block_x, ini_block_y, ini_block_z;
   int i_t, i_block;
   int cpu_number, ng, ratio, n_level, nc;
   int nt, n_block, block_id, lv, n_dim, sys, int_tmp;
   int *p_pos;
   int *n_block_all;
   int *id;
   int *dim;
   int *mesh_type;
   int *level_type, *level_len, *level_index, *level_id, **level_data;
   int i_scalar, n_scalar, *scalar_list, **scalar_type;
   int i_vector, n_vector, *vector_list, **vector_type;

   int var_flag[n_var_total];

   double dx, dy, dz, dx2, dy2, dz2;
   double rangex, rangey, rangez;
   double zero_posx, zero_posy, zero_posz;
   double gm, cfl, gmm1, eps;
   double *t, *dt;
   double *x, *y, *z, *xh, *yh, *zh, *xb, *yb, *zb, **coord;

   #include "data_manage_var_define.include"

   VAR_TYPE *ro, *mx, *my, *mz, *bx, *by, *bz, *en, *gp, *gx, *gy, *gz,
            *roh, *gph, *mph, *teh, *betah, *vxh, *vyh, *vzh, *bxh, *byh, *bzh,
            *jxh, *jyh, *jzh, *fgxh, *fgyh, *fgzh, *fggxh, *fggyh, *fggzh,
            *fmgxh, *fmgyh, *fmgzh, *fmtxh, *fmtyh, *fmtzh;
   VAR_TYPE *rob, *gpb, *mpb, *teb, *betab, *vxb, *vyb, *vzb, *bxb, *byb, *bzb,
            *jxb, *jyb, *jzb, *fgxb, *fgyb, *fgzb, *fggxb, *fggyb, *fggzb,
            *fmgxb, *fmgyb, *fmgzb, *fmtxb, *fmtyb, *fmtzb;
   VAR_TYPE **vector;

   int var_total_dim[n_var_total] = { 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 3, 1, 1, 1, 3, 
                                      1, 1, 1, 3, 1, 1, 1, 3, 1, 1, 1, 3, 1, 1, 1, 3 };
   char var_total_name_list[n_var_total][128] =  { "ro", "gp", "mp", "te", "beta", "vx", "vy", "vz", "vv", 
      "bx", "by", "bz", "bv", "jx", "jy", "jz", "jv", "fgx", "fgy", "fgz", "fgv", "fggx", "fggy", "fggz", "fggv", 
      "fmgx", "fmgy", "fmgz", "fmgv", "fmtx", "fmty", "fmtz", "fmtv" };
   int cpu_id, cpu_n, n_out;


//--------------------------------------------------------------------------------------------------------------------------------//
// start
//--------------------------------------------------------------------------------------------------------------------------------//
   MPI_Init(&argc, &argv);
   MPI_Comm_rank(MPI_COMM_WORLD, &cpu_id);
   MPI_Comm_size(MPI_COMM_WORLD, &cpu_n);

   if (cpu_id == 0)
   {
      printf("\n");
      printf("|------------------|\n");
      printf(" Data manage starts:\n");
      printf("|------------------|\n");
      printf("\n");
   }

   strcpy(char_tmp, "mkdir -p ");
   system(strcat(char_tmp, datavis_path));

//--------------------------------------------------------------------------------------------------------------------------------//
// check no reduplicate variable name
//--------------------------------------------------------------------------------------------------------------------------------//
   for (j = 0; j < n_var_total; j++)
   {
      i_var = 0;
      for (i = 0; i < n_var_total; i++)
         if (strstr(var_total_name_list[j], var_total_name_list[i]) != 0)
         {
            if (i == j) continue;
            i_var++;
         }

      if (i_var > 0)
      {
         printf(" ERROR: reduplicate variable name: %s\n", var_total_name_list[j]);
         for (i = 1; i < n_var_total; i++)
         {
            if (strstr(var_total_name_list[j], var_total_name_list[i]) != 0)
            {
               if (i == j) continue;
               printf(" %s <- %s\n", var_total_name_list[j], var_total_name_list[i]);
            }
         }
         printf("\n");
         return 0;
      }
   }

//--------------------------------------------------------------------------------------------------------------------------------//
// prepare
//--------------------------------------------------------------------------------------------------------------------------------//
   for (i_var = 0; i_var < n_var_total; i_var++) var_flag[i_var] = 0;

   len_str = strlen(var_name_list);
/*   printf("%d\n", len_str);*/
/*   printf("%s\n", var_name_list);*/

   for (i = 0; i < len_str; i++)
      var_name_list[i] = (char) tolower(var_name_list[i]);

   n_var = 0;
   for (i_var = 0; i_var < n_var_total; i_var++)
      if (strstr(var_name_list, var_total_name_list[i_var]) != 0)
      {
         var_flag[i_var] = 1;
         n_var++;
      }

   if (cpu_id == 0)
   {
      printf(" Total number of variables is: %d\n", n_var);
      printf("\n");
      printf(" The following variables will be outputted for visualization:\n");
      printf(" %s\n", var_name_list);
      printf("\n");
   }

   var_name = (char **) malloc(n_var * sizeof(char *));
   for (i_var = 0; i_var < n_var; i_var++) var_name[i_var] = (char *) malloc(char_len_128 * sizeof(char));

   n_scalar = 0;
   n_vector = 0;
   j_var = 0;
   for (i_var = 0; i_var < n_var_total; i_var++)
      if (var_flag[i_var] == 1)
      {
         if (var_total_dim[i_var] == 1) n_scalar++;
         if (var_total_dim[i_var] == 3) n_vector++;
         var_name[j_var] = var_total_name_list[i_var];
         j_var++;
      }

//--------------------------------------------------------------------------------------------------------------------------------//
// read parameters
//--------------------------------------------------------------------------------------------------------------------------------//
   strcpy(char_tmp, datamap_path);
   strcat(char_tmp, "params.mapout");
   datamap_unit = fopen(char_tmp, "r");
   fscanf(datamap_unit, "%s", directory);
   fscanf(datamap_unit, "%s %d", char_tmp, &cpu_number);
   fscanf(datamap_unit, "%s %d", char_tmp, &n_dim);
   fscanf(datamap_unit, "%s %d", char_tmp, &sys);
   fscanf(datamap_unit, "%s %d", char_tmp, &ng);
   fscanf(datamap_unit, "%s %d", char_tmp, &ini_block_x);
   fscanf(datamap_unit, "%s %d", char_tmp, &ini_block_y);
   fscanf(datamap_unit, "%s %d", char_tmp, &ini_block_z);
   fscanf(datamap_unit, "%s %d", char_tmp, &nx_no_bnd);
   fscanf(datamap_unit, "%s %d", char_tmp, &nx_block);
   fscanf(datamap_unit, "%s %d", char_tmp, &nx_total);
   fscanf(datamap_unit, "%s %d", char_tmp, &ny_no_bnd);
   fscanf(datamap_unit, "%s %d", char_tmp, &ny_block);
   fscanf(datamap_unit, "%s %d", char_tmp, &ny_total);
   fscanf(datamap_unit, "%s %d", char_tmp, &nz_no_bnd);
   fscanf(datamap_unit, "%s %d", char_tmp, &nz_block);
   fscanf(datamap_unit, "%s %d", char_tmp, &nz_total);
   fscanf(datamap_unit, "%s %d", char_tmp, &ratio);
   fscanf(datamap_unit, "%s %d", char_tmp, &n_level);
   fscanf(datamap_unit, "%s %lf", char_tmp, &gm);
   fscanf(datamap_unit, "%s %lf", char_tmp, &rangex);
   fscanf(datamap_unit, "%s %lf", char_tmp, &rangey);
   fscanf(datamap_unit, "%s %lf", char_tmp, &rangez);
   fscanf(datamap_unit, "%s %lf", char_tmp, &zero_posx);
   fscanf(datamap_unit, "%s %lf", char_tmp, &zero_posy);
   fscanf(datamap_unit, "%s %lf", char_tmp, &zero_posz);
   fscanf(datamap_unit, "%s %lf", char_tmp, &cfl);
   fclose(datamap_unit);

   if (cpu_id == 0)
   {
      printf(" Work space: %s\n", directory);
      printf("\n");
      printf(" Base x grid points= %d\n", nx_total - ng * 2);
      printf(" Base y grid points= %d\n", ny_total - ng * 2);
      printf(" Base z grid points= %d\n", nz_total - ng * 2);
      printf("\n");
   }

//--------------------------------------------------------------------------------------------------------------------------------//
// read the number of times of outputs
//--------------------------------------------------------------------------------------------------------------------------------//
   strcpy(char_tmp, datamap_path);
   strcat(char_tmp, "nt.mapout");
   datamap_unit = fopen(char_tmp, "rb");
   fread(&int_tmp, sizeof(int), 1, datamap_unit);
   fread(&nt, sizeof(int), 1, datamap_unit);
   fread(&int_tmp, sizeof(int), 1, datamap_unit);
   fclose(datamap_unit);

//--------------------------------------------------------------------------------------------------------------------------------//
// ro, gp, mp, te, beta, vx, vy, vz, vv, bx, by, bz, bv, jx, jy, jz, jv, fgx, fgy, fgz, fgv, fggx, fggy, fggz, fggv, 
// fmgx, fmgy, fmgz, fmgv, fmtx, fmty, fmtz, fmtv
//--------------------------------------------------------------------------------------------------------------------------------//
   nxny = nx_block * ny_block;
   nxnynz = nx_block * ny_block * nz_block;
   mesh_name_block = (char *) malloc(char_len_128 * sizeof(char));
   coord_name = (char **) malloc(n_dim * sizeof(char *));
   coord = (double **) malloc(n_dim * sizeof(double *));
   dim = (int *) malloc(n_dim * sizeof(int));
   vector = (VAR_TYPE **) malloc(n_dim * sizeof(VAR_TYPE *));
   vector_comp_name = (char **) malloc(n_dim * sizeof(char *));
   scalar_name = (char ***) malloc(n_scalar * sizeof(char **));
   vector_name = (char ***) malloc(n_vector * sizeof(char **));
   scalar_type = (int **) malloc(n_scalar * sizeof(int *));
   vector_type = (int **) malloc(n_vector * sizeof(int *));
   scalar_list = (int *) malloc(n_scalar * sizeof(int));
   vector_list = (int *) malloc(n_vector * sizeof(int));

   i_scalar = 0;
   i_vector = 0;
   j_var = 0;
   for (i_var = 0; i_var < n_var_total; i_var++)
      if (var_flag[i_var] == 1)
      {
         if (var_total_dim[i_var] == 1)
         {
            scalar_list[i_scalar] = j_var;
            i_scalar++;
         }
         if (var_total_dim[i_var] == 3)
         {
            vector_list[i_vector] = j_var;
            i_vector++;
         }
         j_var++;
      }

   t = (double *) malloc(nt *  sizeof(double));
   dt = (double *) malloc(nt *  sizeof(double));
   n_block_all = (int *) malloc(nt * sizeof(int));
   level_id = (int *) malloc(n_level * sizeof(int));
   level_index = (int *) malloc(n_level * sizeof(int));
   level_data = (int **) malloc(n_level * sizeof(int *));
   level_len = (int *) malloc(n_level * sizeof(int));
   level_type = (int *) malloc(n_level * sizeof(int));
   p_pos = (int *) malloc(n_dim * sizeof(int));
   id = (int *) malloc(n_dim * n_level * sizeof(int));

   x = (double *) malloc(nx_block * sizeof(double));
   y = (double *) malloc(ny_block * sizeof(double));
   z = (double *) malloc(nz_block * sizeof(double));
   xh = (double *) malloc(nx_block * sizeof(double));
   yh = (double *) malloc(ny_block * sizeof(double));
   zh = (double *) malloc(nz_block * sizeof(double));
   ro = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   mx = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   my = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   mz = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   bx = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   by = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   bz = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   en = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   gp = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   gx = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   gy = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   gz = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   roh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   gph = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   mph = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   teh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   betah = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   vxh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   vyh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   vzh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   bxh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   byh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   bzh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   jxh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   jyh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   jzh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   fgxh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   fgyh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   fgzh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   fggxh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   fggyh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   fggzh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   fmgxh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   fmgyh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   fmgzh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   fmtxh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   fmtyh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));
   fmtzh = (VAR_TYPE *) malloc(nxnynz * sizeof(VAR_TYPE));

   nxnyb = (nx_no_bnd + 1) * (ny_no_bnd + 1);
   nxnynzb = (nx_no_bnd + 1) * (ny_no_bnd + 1) * (ny_no_bnd + 1);
   xb = (double *) malloc((nx_no_bnd + 1) * sizeof(double));
   yb = (double *) malloc((ny_no_bnd + 1) * sizeof(double));
   zb = (double *) malloc((nz_no_bnd + 1) * sizeof(double));
   rob = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   gpb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   mpb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   teb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   betab = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   vxb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   vyb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   vzb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   bxb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   byb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   bzb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   jxb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   jyb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   jzb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   fgxb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   fgyb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   fgzb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   fggxb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   fggyb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   fggzb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   fmgxb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   fmgyb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   fmgzb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   fmtxb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   fmtyb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));
   fmtzb = (VAR_TYPE *) malloc(nxnynzb * sizeof(VAR_TYPE));

//--------------------------------------------------------------------------------------------------------------------------------//
   gmm1 = gm - 1.0;
   eps = 1.0e-12;
   mesh_name_block = strdup("amr_mesh");
   if (output_real8_flag == 4) db_var_type = DB_FLOAT; else db_var_type = DB_DOUBLE;

   if (output_time_start < 0) output_time_start = 0;
   if (output_time_end > nt - 1) output_time_end = nt - 1;

   n_out = 0;

   MPI_Barrier(MPI_COMM_WORLD);

   for (i_t = output_time_start; i_t <= output_time_end; i_t++)
   {
      if (i_t != output_time_start + cpu_id + cpu_n * n_out) continue;
      n_out = n_out + 1;

      sprintf(i_t_char, "%04d", i_t);

//--------------------------------------------------------------------------------------------------------------------------------//
//    read t, dt
//--------------------------------------------------------------------------------------------------------------------------------//
      strcpy(char_tmp, datamap_path);
      strcat(strcat(strcat(char_tmp, "info.t"), i_t_char), ".mapout");
      datamap_unit = fopen(char_tmp, "rb");
      fread(&int_tmp, sizeof(int), 1, datamap_unit);
      fread(&t[i_t], sizeof(double), 1, datamap_unit);
      fread(&int_tmp, sizeof(int), 1, datamap_unit);
      fread(&int_tmp, sizeof(int), 1, datamap_unit);
      fread(&dt[i_t], sizeof(double), 1, datamap_unit);
      fread(&int_tmp, sizeof(int), 1, datamap_unit);
      fclose(datamap_unit);

      strcpy(datavis_name, datavis_path);
      strcat(strcat(datavis_name, "data"), i_t_char);

      printf(" Time = %10.3e\n", t[i_t]);

      strcpy(char_tmp, datavis_name);
      strcat(char_tmp, ".silo");

      dbfile_unit = DBCreate(char_tmp, DB_CLOBBER, DB_LOCAL, NULL, DB_PDB);

      block_id = 0;
      for (i_cpu = 0; i_cpu < cpu_number; i_cpu++)
      {
//--------------------------------------------------------------------------------------------------------------------------------//
//       read n_block
//--------------------------------------------------------------------------------------------------------------------------------//
         sprintf(i_cpu_char, "%04d", i_cpu);

         strcpy(char_tmp, datamap_path);
         strcat(strcat(strcat(strcat(strcat(char_tmp, "data.t"), i_t_char), ".cpu"), i_cpu_char), ".mapout");
         datamap_unit = fopen(char_tmp, "rb");
         fread(&int_tmp, sizeof(int), 1, datamap_unit);
         fread(&n_block_all[i_t], sizeof(int), 1, datamap_unit);
         fread(&int_tmp, sizeof(int), 1, datamap_unit);
         fread(&int_tmp, sizeof(int), 1, datamap_unit);
         fread(level_len, sizeof(int), n_level, datamap_unit);
         fread(&int_tmp, sizeof(int), 1, datamap_unit);
         fread(&int_tmp, sizeof(int), 1, datamap_unit);
         fread(&n_block, sizeof(int), 1, datamap_unit);
         fread(&int_tmp, sizeof(int), 1, datamap_unit);
         if (i_cpu == 0)
         {
            mesh_name = (char **) malloc(n_block_all[i_t] * sizeof(char *));
            mesh_type = (int *) malloc(n_block_all[i_t] * sizeof(int));
            for (i_scalar = 0; i_scalar < n_scalar; i_scalar++)
            {
               scalar_name[i_scalar] = (char **) malloc(n_block_all[i_t] * sizeof(char *));
               scalar_type[i_scalar] = (int *) malloc(n_block_all[i_t] * sizeof(int));
            }
            for (i_vector = 0; i_vector < n_vector; i_vector++)
            {
               vector_name[i_vector] = (char **) malloc(n_block_all[i_t] * sizeof(char *));
               vector_type[i_vector] = (int *) malloc(n_block_all[i_t] * sizeof(int));
            }
            for (i_level = 0; i_level < n_level; i_level++)
            {
               level_id[i_level] = i_level;
               level_data[i_level] = (int *) malloc(level_len[i_level] * sizeof(int));
               level_type[i_level] = DB_BLOCKCENT;
               level_index[i_level] = 0;
            }
         }

         for (i_block = 0; i_block < n_block; i_block++)
         {
//--------------------------------------------------------------------------------------------------------------------------------//
//          read x, ro, mx, my, mz, bx, by, bz, en, gx, gy, gz
//--------------------------------------------------------------------------------------------------------------------------------//
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(id, sizeof(int), n_dim * n_level, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&lv, sizeof(int), 1, datamap_unit);
            fread(&nc, sizeof(int), 1, datamap_unit);
            fread(p_pos, sizeof(int), n_dim, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&dx, sizeof(double), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&dy, sizeof(double), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&dz, sizeof(double), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);

            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(x, sizeof(double), nx_block, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(y, sizeof(double), ny_block, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(z, sizeof(double), nz_block, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(ro, sizeof(VAR_TYPE), nxnynz, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(mx, sizeof(VAR_TYPE), nxnynz, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(my, sizeof(VAR_TYPE), nxnynz, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(mz, sizeof(VAR_TYPE), nxnynz, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(bx, sizeof(VAR_TYPE), nxnynz, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(by, sizeof(VAR_TYPE), nxnynz, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(bz, sizeof(VAR_TYPE), nxnynz, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(en, sizeof(VAR_TYPE), nxnynz, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(gx, sizeof(VAR_TYPE), nxnynz, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(gy, sizeof(VAR_TYPE), nxnynz, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);
            fread(gz, sizeof(VAR_TYPE), nxnynz, datamap_unit);
            fread(&int_tmp, sizeof(int), 1, datamap_unit);

            i_e = nx_block - 1;
            j_e = ny_block - 1;
            k_e = nz_block - 1;

            for (i_grid = 0; i_grid < nxnynz; i_grid++)
               gp[i_grid] = (en[i_grid] - 
               (mx[i_grid] * mx[i_grid] + my[i_grid] * my[i_grid] + mz[i_grid] * mz[i_grid]) / ro[i_grid] / (VAR_TYPE) 2.0 - 
               (bx[i_grid] * bx[i_grid] + by[i_grid] * by[i_grid] + bz[i_grid] * bz[i_grid]) / (VAR_TYPE) 2.0) * (VAR_TYPE) gmm1;

            for (i = 0; i < i_e; i++)
            {
               xh[i] = (x[i] + x[i + 1]) / 2.0;
               if (i >= ng - 1 && i <= nx_block - ng - 1)
                 xb[i - ng + 1] = xh[i];
            }

            for (j = 0; j < j_e; j++)
            {
               yh[j] = (y[j] + y[j + 1]) / 2.0;
               if (j >= ng - 1 && j <= ny_block - ng - 1)
                 yb[j - ng + 1] = yh[j];
            }

            for (k = 0; k < k_e; k++)
            {
               zh[k] = (z[k] + z[k + 1]) / 2.0;
               if (k >= ng - 1 && k <= nz_block - ng - 1)
                 zb[k - ng + 1] = zh[k];
            }

//--------------------------------------------------------------------------------------------------------------------------------//
//          ro, gp, mp, te, beta, vx, vy, vz, vv, bx, by, bz, bv, jx, jy, jz, jv, fgx, fgy, fgz, fgv, fggx, fggy, fggz, fggv, 
//          fmgx, fmgy, fmgz, fmgv, fmtx, fmty, fmtz, fmtv
//--------------------------------------------------------------------------------------------------------------------------------//
            for (i_grid = 0; i_grid < nxnynz; i_grid++)
            {
               k = i_grid / nxny;
               j = (i_grid - k * nxny) / nx_block;
               i = i_grid - k * nxny - j * nx_block;
               if (i >= i_e) continue;
               if (j >= j_e) continue;
               if (k >= k_e) continue;
               p1 = i_grid + 1;
               p2 = i_grid + nx_block;
               p3 = i_grid + nx_block + 1;
               p4 = i_grid + nxny;
               p5 = i_grid + nxny + 1;
               p6 = i_grid + nxny + nx_block;
               p7 = i_grid + nxny + nx_block + 1;
               roh[i_grid] = (ro[i_grid] + ro[p1] + ro[p2] + ro[p3] + ro[p4] + ro[p5] + ro[p6] + ro[p7]) / (VAR_TYPE) 8.0;
               vxh[i_grid] = (mx[i_grid] + mx[p1] + mx[p2] + mx[p3] + mx[p4] + mx[p5] + mx[p6] + mx[p7]) / 
                  (VAR_TYPE) 8.0 / roh[i_grid];
               vyh[i_grid] = (my[i_grid] + my[p1] + my[p2] + my[p3] + my[p4] + my[p5] + my[p6] + my[p7]) / 
                  (VAR_TYPE) 8.0 / roh[i_grid];
               vzh[i_grid] = (mz[i_grid] + mz[p1] + mz[p2] + mz[p3] + mz[p4] + mz[p5] + mz[p6] + mz[p7]) / 
                  (VAR_TYPE) 8.0 / roh[i_grid];
               bxh[i_grid] = (bx[i_grid] + bx[p1] + bx[p2] + bx[p3] + bx[p4] + bx[p5] + bx[p6] + bx[p7]) / (VAR_TYPE) 8.0;
               byh[i_grid] = (by[i_grid] + by[p1] + by[p2] + by[p3] + by[p4] + by[p5] + by[p6] + by[p7]) / (VAR_TYPE) 8.0;
               bzh[i_grid] = (bz[i_grid] + bz[p1] + bz[p2] + bz[p3] + bz[p4] + bz[p5] + bz[p6] + bz[p7]) / (VAR_TYPE) 8.0;
               gph[i_grid] = (gp[i_grid] + gp[p1] + gp[p2] + gp[p3] + gp[p4] + gp[p5] + gp[p6] + gp[p7]) / (VAR_TYPE) 8.0;
               mph[i_grid] = (bxh[i_grid] * bxh[i_grid] + byh[i_grid] * byh[i_grid] + bzh[i_grid] * bzh[i_grid]) / (VAR_TYPE) 2.0;
               teh[i_grid] = gph[i_grid] / roh[i_grid];
               betah[i_grid] = gph[i_grid] / (mph[i_grid] + (VAR_TYPE) eps);
               fgxh[i_grid] = (gx[i_grid] + gx[p1] + gx[p2] + gx[p3] + gx[p4] + gx[p5] + gx[p6] + gx[p7]) / (VAR_TYPE) 8.0;
               fgyh[i_grid] = (gy[i_grid] + gy[p1] + gy[p2] + gy[p3] + gy[p4] + gy[p5] + gy[p6] + gy[p7]) / (VAR_TYPE) 8.0;
               fgzh[i_grid] = (gz[i_grid] + gz[p1] + gz[p2] + gz[p3] + gz[p4] + gz[p5] + gz[p6] + gz[p7]) / (VAR_TYPE) 8.0;
               if (i >= ng - 1 && i <= nx_block - ng - 1)
               if (j >= ng - 1 && j <= ny_block - ng - 1)
               if (k >= ng - 1 && k <= nz_block - ng - 1)
               {
                  i_tmp = i - ng + 1;
                  j_tmp = j - ng + 1;
                  k_tmp = k - ng + 1;
                  i_grid_tmp = k_tmp * nxnyb + j_tmp * (nx_no_bnd + 1) + i_tmp;
                  rob[i_grid_tmp] = roh[i_grid];
                  vxb[i_grid_tmp] = vxh[i_grid];
                  vyb[i_grid_tmp] = vyh[i_grid];
                  vzb[i_grid_tmp] = vzh[i_grid];
                  bxb[i_grid_tmp] = bxh[i_grid];
                  byb[i_grid_tmp] = byh[i_grid];
                  bzb[i_grid_tmp] = bzh[i_grid];
                  gpb[i_grid_tmp] = gph[i_grid];
                  mpb[i_grid_tmp] = mph[i_grid];
                  teb[i_grid_tmp] = teh[i_grid];
                  betab[i_grid_tmp] = betah[i_grid];
                  fgxb[i_grid_tmp] = fgxh[i_grid];
                  fgyb[i_grid_tmp] = fgyh[i_grid];
                  fgzb[i_grid_tmp] = fgzh[i_grid];
               }
            }

            dx2 = dx * 2.0;
            dy2 = dy * 2.0;
            dz2 = dz * 2.0;

            for (i_grid = 0; i_grid < nxnynz; i_grid++)
            {
               k = i_grid / nxny;
               j = (i_grid - k * nxny) / nx_block;
               i = i_grid - k * nxny - j * nx_block;
               if (i == 0 || i >= i_e - 1) continue;
               if (j == 0 || j >= j_e - 1) continue;
               if (k == 0 || k >= k_e - 1) continue;
               p1 = i_grid + 1;
               p2 = i_grid - 1;
               p3 = i_grid + nx_block;
               p4 = i_grid - nx_block;
               p5 = i_grid + nxny;
               p6 = i_grid - nxny;
               jxh[i_grid] = (bzh[p3] - bzh[p4]) / (VAR_TYPE) dy2 - (byh[p5] - byh[p6]) / (VAR_TYPE) dz2;
               jyh[i_grid] = (bxh[p5] - bxh[p6]) / (VAR_TYPE) dz2 - (bzh[p1] - bzh[p2]) / (VAR_TYPE) dx2;
               jzh[i_grid] = (byh[p1] - byh[p2]) / (VAR_TYPE) dx2 - (bxh[p3] - bxh[p4]) / (VAR_TYPE) dy2;
               fggxh[i_grid] = (mph[p1] - mph[p2]) / (VAR_TYPE) dx2;
               fggyh[i_grid] = (mph[p3] - mph[p4]) / (VAR_TYPE) dy2;
               fggzh[i_grid] = (mph[p5] - mph[p6]) / (VAR_TYPE) dz2;
               fmgxh[i_grid] = (mph[p1] - mph[p2]) / (VAR_TYPE) dx2;
               fmgyh[i_grid] = (mph[p3] - mph[p4]) / (VAR_TYPE) dy2;
               fmgzh[i_grid] = (mph[p5] - mph[p6]) / (VAR_TYPE) dz2;
               fmtxh[i_grid] = bxh[i_grid] * (bxh[p1] - bxh[p2]) / (VAR_TYPE) dx2 +
                               byh[i_grid] * (bxh[p3] - bxh[p4]) / (VAR_TYPE) dy2 +
                               bzh[i_grid] * (bxh[p5] - bxh[p6]) / (VAR_TYPE) dz2;
               fmtyh[i_grid] = bxh[i_grid] * (byh[p1] - byh[p2]) / (VAR_TYPE) dx2 +
                               byh[i_grid] * (byh[p3] - byh[p4]) / (VAR_TYPE) dy2 +
                               bzh[i_grid] * (byh[p5] - byh[p6]) / (VAR_TYPE) dz2;
               fmtxh[i_grid] = bxh[i_grid] * (bzh[p1] - bzh[p2]) / (VAR_TYPE) dx2 +
                               byh[i_grid] * (bzh[p3] - bzh[p4]) / (VAR_TYPE) dy2 +
                               bzh[i_grid] * (bzh[p5] - bzh[p6]) / (VAR_TYPE) dz2;
               if (i >= ng - 1 && i <= nx_block - ng - 1)
               if (j >= ng - 1 && j <= ny_block - ng - 1)
               if (k >= ng - 1 && k <= nz_block - ng - 1)
               {
                  i_tmp = i - ng + 1;
                  j_tmp = j - ng + 1;
                  k_tmp = k - ng + 1;
                  i_grid_tmp = k_tmp * nxnyb + j_tmp * nx_no_bnd + i_tmp;
                  jxb[i_grid_tmp] = jxh[i_grid];
                  jyb[i_grid_tmp] = jyh[i_grid];
                  jzb[i_grid_tmp] = jzh[i_grid];
                  fggxb[i_grid_tmp] = fggxh[i_grid];
                  fggyb[i_grid_tmp] = fggyh[i_grid];
                  fggzb[i_grid_tmp] = fggzh[i_grid];
                  fmgxb[i_grid_tmp] = fmgxh[i_grid];
                  fmgyb[i_grid_tmp] = fmgyh[i_grid];
                  fmgzb[i_grid_tmp] = fmgzh[i_grid];
                  fmtxb[i_grid_tmp] = fmtxh[i_grid];
                  fmtyb[i_grid_tmp] = fmtyh[i_grid];
                  fmtzb[i_grid_tmp] = fmtzh[i_grid];
               }
            }

//--------------------------------------------------------------------------------------------------------------------------------//
//          silo ouput
//--------------------------------------------------------------------------------------------------------------------------------//
            sprintf(block_id_char, "/block%d", block_id);
            mesh_name[block_id] = strdup("");
            strcpy(mesh_name[block_id], block_id_char);
            strcat(strcat(mesh_name[block_id], "/"), mesh_name_block);
            mesh_type[block_id] = DB_QUADMESH;
/*            printf("%s\n", block_id_char);*/
            for (i_scalar = 0; i_scalar < n_scalar; i_scalar++)
            {
               scalar_name[i_scalar][block_id] = strdup("");
               strcpy(scalar_name[i_scalar][block_id], block_id_char);
               strcat(strcat(scalar_name[i_scalar][block_id], "/"), var_name[scalar_list[i_scalar]]);
               scalar_type[i_scalar][block_id] = DB_QUADVAR;
            }
            for (i_vector = 0; i_vector < n_vector; i_vector++)
            {
               vector_name[i_vector][block_id] = strdup("");
               strcpy(vector_name[i_vector][block_id], block_id_char);
               strcat(strcat(vector_name[i_vector][block_id], "/"), var_name[vector_list[i_vector]]);
               vector_type[i_vector][block_id] = DB_QUADVAR;
            }
            level_data[lv - 1][level_index[lv - 1]++] = block_id;
            DBMkDir(dbfile_unit, block_id_char);
            DBSetDir(dbfile_unit, block_id_char);
            coord_name[0] = strdup("x");
            coord_name[1] = strdup("y");
            coord_name[2] = strdup("z");
            coord[0] = xb;
            coord[1] = yb;
            coord[2] = zb;
            dim[0] = nx_no_bnd + 1;
            dim[1] = ny_no_bnd + 1;
            dim[2] = nz_no_bnd + 1;
            DBPutQuadmesh(dbfile_unit, mesh_name_block, coord_name, coord, dim, n_dim, DB_DOUBLE, DB_COLLINEAR, NULL);
            i_var = 0;
            j_var = 0;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, rob, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, gpb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, mpb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, teb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, betab, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, vxb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, vyb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, vzb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               vector_comp_name[0] = strdup("vx");
               vector_comp_name[1] = strdup("vy");
               vector_comp_name[2] = strdup("vz");
               vector[0] = vxb;
               vector[1] = vyb;
               vector[2] = vzb;
               DBPutQuadvar(dbfile_unit, var_name[j_var], mesh_name_block, n_dim, vector_comp_name, vector, dim, n_dim, NULL, 0,
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, bxb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, byb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, bzb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               vector_comp_name[0] = strdup("bx");
               vector_comp_name[1] = strdup("by");
               vector_comp_name[2] = strdup("bz");
               vector[0] = bxb;
               vector[1] = byb;
               vector[2] = bzb;
               DBPutQuadvar(dbfile_unit, var_name[j_var], mesh_name_block, n_dim, vector_comp_name, vector, dim, n_dim, NULL, 0,
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, jxb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, jyb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, jzb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               vector_comp_name[0] = strdup("jx");
               vector_comp_name[1] = strdup("jy");
               vector_comp_name[2] = strdup("jz");
               vector[0] = jxb;
               vector[1] = jyb;
               vector[2] = jzb;
               DBPutQuadvar(dbfile_unit, var_name[j_var], mesh_name_block, n_dim, vector_comp_name, vector, dim, n_dim, NULL, 0,
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, fgxb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, fgyb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, fgzb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               vector_comp_name[0] = strdup("fgx");
               vector_comp_name[1] = strdup("fgy");
               vector_comp_name[2] = strdup("fgz");
               vector[0] = fgxb;
               vector[1] = fgyb;
               vector[2] = fgzb;
               DBPutQuadvar(dbfile_unit, var_name[j_var], mesh_name_block, n_dim, vector_comp_name, vector, dim, n_dim, NULL, 0,
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, fggxb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, fggyb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, fggzb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               vector_comp_name[0] = strdup("fggx");
               vector_comp_name[1] = strdup("fggy");
               vector_comp_name[2] = strdup("fggz");
               vector[0] = fggxb;
               vector[1] = fggyb;
               vector[2] = fggzb;
               DBPutQuadvar(dbfile_unit, var_name[j_var], mesh_name_block, n_dim, vector_comp_name, vector, dim, n_dim, NULL, 0,
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, fmgxb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, fmgyb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, fmgzb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               vector_comp_name[0] = strdup("fmgx");
               vector_comp_name[1] = strdup("fmgy");
               vector_comp_name[2] = strdup("fmgz");
               vector[0] = fmgxb;
               vector[1] = fmgyb;
               vector[2] = fmgzb;
               DBPutQuadvar(dbfile_unit, var_name[j_var], mesh_name_block, n_dim, vector_comp_name, vector, dim, n_dim, NULL, 0,
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, fmtxb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, fmtyb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               DBPutQuadvar1(dbfile_unit, var_name[j_var], mesh_name_block, fmtzb, dim, n_dim, NULL, 0, 
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            i_var++;
            if (var_flag[i_var] == 1)
            {
               vector_comp_name[0] = strdup("fmtx");
               vector_comp_name[1] = strdup("fmty");
               vector_comp_name[2] = strdup("fmtz");
               vector[0] = fmtxb;
               vector[1] = fmtyb;
               vector[2] = fmtzb;
               DBPutQuadvar(dbfile_unit, var_name[j_var], mesh_name_block, n_dim, vector_comp_name, vector, dim, n_dim, NULL, 0,
                  db_var_type, DB_NODECENT, NULL);
               j_var++;
            }
            DBSetDir(dbfile_unit, "..");
            block_id++;
         }
      }
      opt_list = DBMakeOptlist(10);
      DBAddOption(opt_list, DBOPT_CYCLE, &i_t);
      DBAddOption(opt_list, DBOPT_DTIME, &t[i_t]);
      DBAddOption(opt_list, DBOPT_MRGTREE_NAME, "mrgtree");
      DBPutMultimesh(dbfile_unit, mesh_name_block, n_block_all[i_t], mesh_name, mesh_type, opt_list);
/*for (i_block = 0; i_block < n_block_all[i_t]; i_block++) printf("%s %d\n", mesh_name[i_block], mesh_type[i_block]);*/
      for (i_scalar = 0; i_scalar < n_scalar; i_scalar++)
         DBPutMultivar(dbfile_unit, var_name[scalar_list[i_scalar]], n_block_all[i_t], 
            scalar_name[i_scalar], scalar_type[i_scalar], 0);
      for (i_vector = 0; i_vector < n_vector; i_vector++)
         DBPutMultivar(dbfile_unit, var_name[vector_list[i_vector]], n_block_all[i_t], 
            vector_name[i_vector], vector_type[i_vector], 0);

      DBPutGroupelmap(dbfile_unit, "amr_levels", n_level, level_type, level_len, 0, level_data, 0, 0, 0);
      mrgtree = DBMakeMrgtree(DB_MULTIMESH, 0, 1, 0);
      DBAddRegion(mrgtree, "amr_decomp", 0, 2, 0, 0, 0, 0, 0, 0); 
      DBSetCwr(mrgtree, "amr_decomp");
      DBAddRegion(mrgtree, "levels", 0, n_level, 0, 0, 0, 0, 0, 0);
      DBSetCwr(mrgtree, "levels");
      level_name[0] = "@level%d@n";
      DBAddRegionArray(mrgtree, n_level, level_name, 0, "amr_levels", 1, level_id, level_len, level_type, 0);

      DBSetCwr(mrgtree, "..");
      DBAddRegion(mrgtree, "patches", 0, n_block_all[i_t], 0, 0, 0, 0, 0, 0);
      DBSetCwr(mrgtree, "..");
      DBSetCwr(mrgtree, "..");
      DBPutMrgtree(dbfile_unit, "mrgtree", mesh_name_block, mrgtree, 0);

      DBFreeMrgtree(mrgtree);
      DBFreeOptlist(opt_list);
      DBClose(dbfile_unit);
      free(mesh_name);
      free(mesh_type);
      for (i_scalar = 0; i_scalar < n_scalar; i_scalar++)
      {
         free(scalar_name[i_scalar]);
         free(scalar_type[i_scalar]);
      }
      for (i_vector = 0; i_vector < n_vector; i_vector++)
      {
         free(vector_name[i_vector]);
         free(vector_type[i_vector]);
      }
      for (i_level = 0; i_level < n_level; i_level++)
         free(level_data[i_level]);
   }

//--------------------------------------------------------------------------------------------------------------------------------//
// free memory
//--------------------------------------------------------------------------------------------------------------------------------//
   free(t);
   free(dt);
   free(n_block_all);
   free(level_id);
   free(level_index);
   free(level_data);
   free(level_len);
   free(level_type);
   free(p_pos);
   free(id);

   free(x);
   free(y);
   free(z);
   free(xh);
   free(yh);
   free(zh);
   free(ro);
   free(mx);
   free(my);
   free(mz);
   free(bx);
   free(by);
   free(bz);
   free(en);
   free(gp);
   free(gx);
   free(gy);
   free(gz);
   free(roh);
   free(gph);
   free(mph);
   free(teh);
   free(betah);
   free(vxh);
   free(vyh);
   free(vzh);
   free(bxh);
   free(byh);
   free(bzh);
   free(jxh);
   free(jyh);
   free(jzh);
   free(fgxh);
   free(fgyh);
   free(fgzh);
   free(fggxh);
   free(fggyh);
   free(fggzh);
   free(fmgxh);
   free(fmgyh);
   free(fmgzh);
   free(fmtxh);
   free(fmtyh);
   free(fmtzh);

   free(xb);
   free(yb);
   free(zb);
   free(rob);
   free(gpb);
   free(mpb);
   free(teb);
   free(betab);
   free(vxb);
   free(vyb);
   free(vzb);
   free(bxb);
   free(byb);
   free(bzb);
   free(jxb);
   free(jyb);
   free(jzb);
   free(fgxb);
   free(fgyb);
   free(fgzb);
   free(fggxb);
   free(fggyb);
   free(fggzb);
   free(fmgxb);
   free(fmgyb);
   free(fmgzb);
   free(fmtxb);
   free(fmtyb);
   free(fmtzb);

   MPI_Barrier(MPI_COMM_WORLD);

   if (cpu_id == 0)
   {
      printf("\n");
      printf("|--------------------|\n");
      printf(" Data manage finished.\n");
      printf("|--------------------|\n");
      printf("\n");
   }

   MPI_Finalize();
//--------------------------------------------------------------------------------------------------------------------------------//
   return 0;
}
